Data processing

Process the capture data

The raw capture data is loaded and processed as described in the manuscript. The basic processing steps are:
- retain all first capture records (these are all pre-treatment)
- for second and third captures, remove any that occur after manipulative experiments (cort dosing, predators, taping)
- subset to captures where we actually have at least one cort measure
- some cleanup for erroneous values and data wrangling (not changing any data or subsetting here) - in the end, this results in the inclusion of 1428 rows for baseline, 1307 rows for stress-induced, and 1238 rows for dex

## Load data ----
  # Load data for each birds capture
      db <- read.delim(here::here("1_raw_data/data_by_capture.txt"))
      # Subset to non invasive experiments
        db <- subset(db, db$ok_for_weather == "yes")
      # subset to females
        db <- subset(db, db$sex == "F")
        
  # add in capture hour column
      # Capture hour
        db$cap_hr <- round(db$cap_time / 60, 0)
        
  # fix some cort values being in different columns and correct some values for extraction efficiency
            for(i in 1:nrow(db)){
              if(is.na(db$cort1_correct[i]) == TRUE){
                if(is.na(db$cort1_raw[i]) == FALSE){db$cort1_correct[i] <- db$cort1_raw[i] * (100/85)}
              }
              if(is.na(db$cort2_correct[i]) == TRUE){
                if(is.na(db$cort1_raw[i]) == FALSE){db$cort2_correct[i] <- db$cort2_raw[i] * (100/85)}
              }
              if(is.na(db$cort3_correct[i]) == TRUE){
                if(is.na(db$cort1_raw[i]) == FALSE){db$cort3_correct[i] <- db$cort3_raw[i] * (100/85)}
              }
            }
        
  # subset to captures where we have at least one cort measurement      
        db <- subset(db, db$cort1_correct > 0 | db$cort2_correct > 0 | db$cort3_correct > 0)
        db_save <- db # saving unmodified version here

Process the weather data

I’m using hourly weather data from the game farm road weather station accessed with from the Northeast Regional Climate Center. The file includes records from January 1983 to July 2020. Some basic filtering is done to remove very high or low values for temperature or time (these are likely either error codes like ‘999’ or actual data errors). I did a bunch of descriptive plots just to look for mistakes in the data. This represents 38 years of weather data. The raw data should be accessed from the Northeast Climate Center directly, but the repository includes the code here to illustrate the processing along with the capture data after summarized weather data has been merged.

  # Load weather data
      dw <- read.delim(here::here("1_raw_data/Input_Weather_Data.txt"))      ## NEED TO ACCESS DATA FROM NORTHEAST CLIMATE CENTER
    
      # Subset missing data values or errors
        dw <- subset(dw, dw$avg_temp_C < 200 & dw$hour < 25 & dw$avg_temp_C > -100)
        
  # filter to daytime
        dw_day <- dw %>% filter(hour > 5, hour < 21)
        dw_day <- dw_day %>% filter(doy > 90 & doy < 215)
        
  # this has the daytime average temperature for every year day combo to be used later for calculating expected temperature      
        dw_day_mu <- dw_day %>%
          group_by(year, doy) %>%
          dplyr::summarise(av_temp_C = mean(avg_temp_C, na.rm = TRUE), n_records = n(), hi_temp_C = max(avg_temp_C, na.rm = TRUE))
        dw_day_mu$avg_temp_C <- dw_day_mu$av_temp_C

Merge capture and weather data

Joining the capture records to weather data. First I’m matching weather for the sliding window analysis using only the actual weather (not long term expectations). The code is below, but the basic process is:

  • loop through each capture record

  • subset out temperature records from daytime (6am to 8pm)

  • subset to the date range; date ranges start 1-30 days before capture and extend from 1 day long up to the day of capture. For starting 30 days before capture there are 30 date ranges that are 1, 2, 3…29, 30 days long.

  • calculate the mean temperature for the subset window and add it to the capture data frame

  • We decided to limit this range to 7 days for base and stress and 21 days for dex, since I already had the code written I’m still calculating all of them, but I’ve clipped it down to just those ranges.

Next, expected temperatures in a separate loop, but only for the time periods we identified in the sliding window analysis

# ## Add weather data to each capture
#   # Creating a list of labels for each time period that weather will be calculated for
    lab_cols <- paste("tb", seq(1, 30, 1), sep = "_")
    for(i in 1:30){
      lab_cols2 <- paste(lab_cols[i:30], rep(i, 31 - i), sep = "_")
      if(i == 1){lab_cols3 <- lab_cols2}
      if(i > 1){lab_cols3 <- c(lab_cols3, lab_cols2)}
    }
# 
#     dblist <- as.list(NA)
#     dblist_ex <- as.list(NA)
#     dblist_sd <- as.list(NA)
# 
#   # Looping through every capture record
#     for(i in 1:nrow(db)){
#       # Creating subsets for different time periods that were predefined before we switched to sliding window
#             sub <- subset(dw, dw$year == db$year[i])
#             sub1 <- subset(sub, sub$doy == db$cap_doy[i] & sub$hour < db$cap_hr[i] + 1 & sub$hour > db$cap_hr[i] - 3)
#             sub1b <- subset(sub, sub$doy == db$cap_doy[i] & sub$hour < 9 & sub$hour > 5)
#             sub2a <- subset(sub, sub$doy == db$cap_doy[i] | sub$doy == db$cap_doy[i] - 1)
#             sub2 <- subset(sub2a, sub2a$hour > 21 | sub2a$hour < 6)
#             sub3 <- subset(sub, sub$doy == db$cap_doy[i] - 1 & sub$hour > 5 & sub$hour < 21)
#             sub4 <- subset(sub, sub$doy < db$cap_doy[i] & sub$doy > db$cap_doy[i] - 4 & sub$hour > 5 & sub$hour < 21)
#             sub5 <- subset(sub, sub$doy < db$cap_doy[i] & sub$doy > db$cap_doy[i] - 8 & sub$hour > 5 & sub$hour < 21)
#             sub6 <- subset(sub, sub$doy < db$cap_doy[i] - 3 & sub$doy > db$cap_doy[i] - 7 & sub$hour > 5 & sub$hour < 21)
#             sub7 <- subset(sub, sub$doy < db$cap_doy[i] - 6 & sub$doy > db$cap_doy[i] - 10 & sub$hour > 5 & sub$hour < 21)
#             sub8 <- subset(sub, sub$doy < db$cap_doy[i] - 9 & sub$doy > db$cap_doy[i] - 13 & sub$hour > 5 & sub$hour < 21)
#             sub9 <- subset(sub, sub$doy < db$cap_doy[i] - 12 & sub$doy > db$cap_doy[i] - 26 & sub$hour > 5 & sub$hour < 21)
#             sub10 <- subset(sub, sub$doy < db$cap_doy[i] - 7 & sub$doy > db$cap_doy[i] - 9 & sub$hour > 5 & sub$hour < 21)
#             sub11 <- subset(sub, sub$doy == db$cap_doy[i] & sub$hour > 5 & sub$hour < 21)
# 
#       # Adding the predefined temperature records
#             db$t_3cap_av[i] <- mean(na.omit(sub1$avg_temp_C))
#             db$t_m_av[i] <- mean(na.omit(sub1b$avg_temp_C))
#             db$t_nb_av[i] <- mean(na.omit(sub2$avg_temp_C))
#             db$t_nb_lo[i] <- min(na.omit(sub2$avg_temp_C))
#             db$t_db_av[i] <- mean(na.omit(sub3$avg_temp_C))
#             db$t_db_hi[i] <- max(na.omit(sub3$avg_temp_C))
#             db$t_3b_av[i] <- mean(na.omit(sub4$avg_temp_C))
#             db$t_3b_hi[i] <- max(na.omit(sub4$avg_temp_C))
#             db$t_7b_av[i] <- mean(na.omit(sub5$avg_temp_C))
#             db$t_7b_hi[i] <- max(na.omit(sub5$avg_temp_C))
#             db$t_4to6_av[i] <- mean(na.omit(sub6$avg_temp_C))
#             db$t_7to9_av[i] <- mean(na.omit(sub7$avg_temp_C))
#             db$t_10to12_av[i] <- mean(na.omit(sub8$avg_temp_C))
#             db$t_13to25_av[i] <- mean(na.omit(sub9$avg_temp_C))
#             db$t_7_av[i] <- mean(na.omit(sub10$avg_temp_C))
#             db$min_cap_day[i] <- min(na.omit(sub11$avg_temp_C))
# 
#       # Nested loop that is calculating the sliding window temperature records. These are stored in their own object.
#           # Right now it records daytime temperature for each day starting 1-30 days ago and with the averaging period
#             # extending from 1 day up through the day of capture.
#             subx <- subset(sub, sub$hour > 5 & sub$hour < 21)
#             subxs <- subx %>%
#               group_by(doy) %>%
#               dplyr::summarise(avg_temp_C = mean(avg_temp_C, na.rm = TRUE), max_temp_C = max(avg_temp_C, na.rm = TRUE))
#             ss1 <- subxs[which(subxs$doy == db$cap_doy[i] - 1):which(subxs$doy == db$cap_doy[i] - 30), "avg_temp_C"]
#             ss2 <- rep(NA, 465)
#             ss2[1:30] <- ss1$avg_temp_C
#             counter <- 31
#             for(t in 1:29){
#               if(t == 1){
#                 ss1x <- data.frame(V1 = ss1$avg_temp_C, V2 = c(ss1$avg_temp_C[(1+t):30], rep(NA, t)))
#               }
#               if(t > 1){
#                 ss1x[, t+1] <- c(ss1$avg_temp_C[(1+t):30], rep(NA, t))
#               }
#             }
#             for(t in 1:29){
#               ss1x2 <- ss1x[1:(30-t), 1:(t + 1)]
#               ss2[counter:(counter + 29 - t)] <- rowMeans(ss1x2, na.rm = TRUE)
#               counter <- counter + 30 - t
#             }
# 
#             dblist[[i]] <- ss2
# 
#             print(paste0(i, " of ", nrow(db)))
#             
#         # Nested loop similar to above but calculating offset from expected NEEDS WORK
#             dw_day_xs <- dw_day_mu %>%
#               group_by(doy) %>%
#               dplyr::summarise(avg_temp_C2 = mean(avg_temp_C, na.rm = TRUE), avg_temp_sd = sd(avg_temp_C, na.rm = TRUE), n = n())
#             ss1_ex <- dw_day_xs[which(dw_day_xs$doy == db$cap_doy[i] - 1):which(dw_day_xs$doy == db$cap_doy[i] - 30), "avg_temp_C2"]
#             ss1_sd <- dw_day_xs[which(dw_day_xs$doy == db$cap_doy[i] - 1):which(dw_day_xs$doy == db$cap_doy[i] - 30), "avg_temp_sd"]
#             ss2_ex <- rep(NA, 465)
#             ss2_sd <- rep(NA, 465)
#             ss2_ex[1:30] <- ss1_ex$avg_temp_C2
#             ss2_sd[1:30] <- ss1_sd$avg_temp_sd
#           
#         # Loop for expected average temperature for the date range      
#             counter <- 31
#             for(t in 1:29){
#               if(t == 1){
#                 ss1x_ex <- data.frame(V1 = ss1_ex$avg_temp_C2, V2 = c(ss1_ex$avg_temp_C2[(1+t):30], rep(NA, t)))
#               }
#               if(t > 1){
#                 ss1x_ex[, t+1] <- c(ss1_ex$avg_temp_C2[(1+t):30], rep(NA, t))
#               }
#             }
#             for(t in 1:29){
#               ss1x2_ex <- ss1x_ex[1:(30-t), 1:(t + 1)]
#               ss2_ex[counter:(counter + 29 - t)] <- rowMeans(ss1x2_ex, na.rm = TRUE)
#               counter <- counter + 30 - t
#             }
# 
#             dblist_ex[[i]] <- ss2_ex
#             
#         # Loop for expected standard deviation for the date range
#             counter <- 31
#             for(t in 1:29){
#               if(t == 1){
#                 ss1x_sd <- data.frame(V1 = ss1_sd$avg_temp_sd, V2 = c(ss1_sd$avg_temp_sd[(1+t):30], rep(NA, t)))
#               }
#               if(t > 1){
#                 ss1x_sd[, t+1] <- c(ss1_sd$avg_temp_sd[(1+t):30], rep(NA, t))
#               }
#             }
#             for(t in 1:29){
#               ss1x2_sd <- ss1x_sd[1:(30-t), 1:(t + 1)]
#               ss2_sd[counter:(counter + 29 - t)] <- rowMeans(ss1x2_sd, na.rm = TRUE)
#               counter <- counter + 30 - t
#             }
# 
#             dblist_sd[[i]] <- ss2_sd
#             
#     }
# 
#     # adding the object with the sliding window temperatures back to the breeding records
#         db[, lab_cols3] <- NA
#         db_dev <- db
#         for(i in 1:nrow(db)){
#           db[i, 52:516] <- dblist[[i]]
#         }
#     
#     # adding the deviations instead, keeping this as a separate object for now
#         db_dev[, lab_cols3] <- NA
        # db_deviation <- as.list(NA)
        # for(i in 1:length(dblist)){
        #   db_deviation[[i]] <- (dblist[[i]] - dblist_ex[[i]])
        #        #/ dblist_sd[[i]]
        # }
        # for(i in 1:nrow(db_dev)){
        #   db_dev[i, 52:516] <- db_deviation[[i]]
        # }
#           

    # Add in expected deviation for day of capture (separate from the other loop above)
      #   sub1 <- subset(dw, dw$hour > 5 & dw$hour < 9)
      #   sub2 <- sub1 %>%
      #     group_by(year, doy) %>%
      #     dplyr::summarise(t_m_av = mean(avg_temp_C))
      # for(u in 1:nrow(db_dev)){
      #   sub3 <- subset(sub2, sub2$doy == db_dev$cap_doy[u])
      #   db_dev$t_m_av[u] <- (db_dev$t_m_av[u] - mean(sub3$t_m_av)) / sd(sub3$t_m_av)
      # }

# 
#     ## Fix some inf values
#       db$t_nb_lo[which(!is.finite(db$t_nb_lo))] <- NA
#       db_dev$t_nb_lo[which(!is.finite(db_dev$t_nb_lo))] <- NA

    # the loop above takes a while to run so I've saved the output and reloaded it here. Should be run again if any changes are needed.
      #saveRDS(db, file = here::here("2_modified_data", "weather_cort.rds"))
      #saveRDS(db_dev, file = here::here("2_modified_data", "weather_cort_ex.rds"))
      db <- readRDS(here::here("2_modified_data", "weather_cort.rds"))
      db_dev <- readRDS(here::here("2_modified_data", "weather_cort_ex.rds"))
      
    # ## Log transform 
       db$cort1_correct <- log(db$cort1_correct)
       db$cort2_correct <- log(db$cort2_correct)
       db$cort3_correct <- log(db$cort3_correct)
       db_dev$cort1_correct <- log(db_dev$cort1_correct)
       db_dev$cort2_correct <- log(db_dev$cort2_correct)
       db_dev$cort3_correct <- log(db_dev$cort3_correct)
      
    ## Separate data for base, stress, and dex
      base_af <- subset(db, db$type1 == "B" & is.na(db$type1) == FALSE)
      stress_af <- subset(db, db$type2 == "S" & is.na(db$type2) == FALSE & is.na(db$type1) == FALSE)
      dex_af <- subset(db, db$type3 == "D" & is.na(db$type3) == FALSE & is.na(db$type2) == FALSE)
      
    ## Separate data for base, stress, and dex
      base_af_ex <- subset(db_dev, db_dev$type1 == "B" & is.na(db_dev$type1) == FALSE)
      stress_af_ex <- subset(db_dev, db_dev$type2 == "S" & is.na(db_dev$type2) == FALSE & is.na(db_dev$type1) == FALSE)
      dex_af_ex <- subset(db_dev, db_dev$type3 == "D" & is.na(db_dev$type3) == FALSE & is.na(db_dev$type2) == FALSE)  

Descriptive plots, sample sizes, etc.

Here I’m just making a few plots and summaries that might be useful for describing the dataset before doing actual analyses.

p1 <- ggplot(data = dw_day_mu, mapping = aes(x = doy, y = avg_temp_C)) + 
      geom_point(alpha = 0.2, color = "coral3") + 
      geom_smooth(fill = "slateblue", color = "slateblue") +
      theme_bw() +
      theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +
      theme(axis.text = element_text(size = 12), axis.title = element_text(size = 16)) +
      ylab("Average daytime temperature \u00B0C") +
      xlab("Day of year") +
      coord_cartesian(xlim = c(135, 185), ylim = c(4, 32)) +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "B", size = 6)
      
p2 <- ggplot(data = db, mapping = aes(x = cap_doy)) +
      geom_histogram(binwidth = 1, fill = "gray25", alpha = 0.4, color = "black") +
      theme_bw() +
      theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +
      theme(axis.text = element_text(size = 12), axis.title = element_text(size = 16)) +
      ylab("Number of captures") +
      xlab("") +
      coord_cartesian(xlim = c(135, 185)) +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "A", size = 6)
  
ggpubr::ggarrange(p2, p1, nrow = 2, heights = c(1, 1.5))
Distribution of capture dates and distribution of average daytime temperatures on each day of the year. Red points are individual days and blue line is the average over 38 years.

Distribution of capture dates and distribution of average daytime temperatures on each day of the year. Red points are individual days and blue line is the average over 38 years.

After the filtering described above, there are a total of 1429 captures included from a total of 719 unique females. There are 892 unique female-year combinations (females caught in different years). Reporting on just the three sample types, there are 1428 separate baseline measurements, 1307 separate stress measurements (these all have matched base measures), and 1238 separate dex measurements (these all have paired stress measures). That means there are a total of 3973 separate corticosterone measurements included. For models fit later on I will need to restrict to records with no missing data (in fixed or random effects) so the sample sizes might differ a bit from the full dataset here.

Here is the distribution by year:

ggplot(db, mapping = aes(x = cap_doy)) + 
  geom_histogram(binwidth = 1) + 
  facet_wrap(~as.factor(year)) +
  theme_bw() +
  xlab("Capture day of year") +
  ylab("Number of captures")

Correlation between date and temperature

This is calculating the correlation between day of year and daytime average temperature. Only days where a capture happened are included and if multiple captures happened on a given day the datapoint is only counted once. The simple pairwise correlation and p-value are printed below.

doycor <- base_af %>%
  dplyr::group_by(year, cap_doy) %>%
  dplyr::summarise(n = n(), temperature = mean(t_db_av))

#cor(doycor$cap_doy, doycor$temperature)
cor.test(doycor$cap_doy, doycor$temperature)
## 
##  Pearson's product-moment correlation
## 
## data:  doycor$cap_doy and doycor$temperature
## t = 5.0823, df = 261, p-value = 7.117e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1858762 0.4062953
## sample estimates:
##       cor 
## 0.3000861

Analyses

1. Does ambient temperature causally affect corticosterone levels

This is from the cold snap experiment in 2021 where nests were either cooled with ice packs or controls.

Description of sampling: females had ice packs in their nest or control packs starting on morning of the 4th day after hatching. Ice packs were swapped at 6/9/12/3/6 each day. Females were captured early in morning on day 7 after hatching and three blood samples were taken at randomly selected time points close to 0/5/10/20/30/60 minutes.

Models are fit as a GAMM and tables here are for the GAMM and the linear part of the GAMM with the random effect.

Figures are produced using the raw data and the GAMM fit.

d_exp <- read.delim(here::here("1_raw_data/cold_snap_adult.txt"))

d_exp1 <- subset(d_exp, d_exp$cap_num == 1)
d_exp1$type <- as.factor(d_exp1$type)
d_exp1$type <- relevel(d_exp1$type, ref = "S")
d_exp1$type <- relevel(d_exp1$type, ref = "B")
p_exp1 <- ggplot(data = d_exp1, mapping = aes(x = type, y = final_cort, color = treatment, fill = treatment)) +
  geom_boxplot(alpha = 0.4) +
  #geom_point() +
  scale_color_manual(values = c("slateblue", "orange")) +
  scale_fill_manual(values = c("slateblue", "orange")) +
  xlab("Sample type") + 
  ylab("Corticosterone (ng/ml)") +
  theme_bw() +
  theme(legend.position = c(0.8, 0.8), axis.title = element_text(size = 12)) +
  coord_cartesian(ylim = c(0, 110))


d_exp2 <- subset(d_exp, d_exp$cap_num == 2)
d_exp2$treatment2 <- "cooled"
for(i in 1:nrow(d_exp2)){
  if(d_exp2$treatment[i] == "Cold_Control"){d_exp2$treatment2[i] <- "control"}
}
d_exp2$treatment2 <- factor(d_exp2$treatment2, levels = c("cooled", "control"))
p_exp2 <- ggplot(data = d_exp2, mapping = aes(x = latency / 60, y = final_cort, color = treatment2, fill = treatment2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "gam") +
  scale_color_manual(values = c("slateblue", "orange")) +
  scale_fill_manual(values = c("slateblue", "orange")) +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
        legend.title = element_blank(), legend.text = element_text(size = 14),
        axis.title = element_text(size = 14), axis.text = element_text(size = 12)) +
  annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "B", size = 6) +
  xlab("Sampling latency (minutes)") +
  ylab("Corticosterone (ng/ml)") +
  theme(legend.position = c(0.75, 0.8)) +
  coord_cartesian(ylim = c(0, 110))

#ggpubr::ggarrange(p1, p2)

d_exp2$lat_min <- d_exp2$latency/60
m_cold <- lmer(final_cort ~ lat_min + I(lat_min^2) + I(lat_min^3) + treatment + (1|band), data = d_exp2)
d_exp2$trt3 <- factor(d_exp2$treatment2, ordered = FALSE)
m_cold2 <- gam(final_cort ~ trt3 + s(lat_min, by = trt3, k = -1), data = d_exp2)
m_cold2b <- gam(final_cort ~ trt3 + s(lat_min, by = trt3, k = -1) + s(band, bs = "re"), data = d_exp2)
m_cold2c <- gamm(final_cort ~ trt3 + s(lat_min, by = trt3, k = -1), data = d_exp2, random = list(band = ~1))

tab_model(m_cold2c$gam)
  final cort
Predictors Estimates CI p
(Intercept) 30.47 26.41 – 34.54 <0.001
trt3control -10.51 -16.61 – -4.41 0.001
Smooth term (lat_min) *
trt3cooled
<0.001
Smooth term (lat_min) *
trt3control
<0.001
Observations 98
R2 0.582
summary(m_cold2c$lme)
## Linear mixed-effects model fit by maximum likelihood
##  Data: strip.offset(mf) 
##        AIC      BIC    logLik
##   790.1552 810.8349 -387.0776
## 
## Random effects:
##  Formula: ~Xr - 1 | g
##  Structure: pdIdnot
##              Xr1      Xr2      Xr3      Xr4      Xr5      Xr6      Xr7      Xr8
## StdDev: 37.87277 37.87277 37.87277 37.87277 37.87277 37.87277 37.87277 37.87277
## 
##  Formula: ~Xr.0 - 1 | g.0 %in% g
##  Structure: pdIdnot
##            Xr.01    Xr.02    Xr.03    Xr.04    Xr.05    Xr.06    Xr.07    Xr.08
## StdDev: 18.65853 18.65853 18.65853 18.65853 18.65853 18.65853 18.65853 18.65853
## 
##  Formula: ~1 | band %in% g.0 %in% g
##         (Intercept) Residual
## StdDev:    6.624314 9.681687
## 
## Fixed effects: y ~ X - 1 
##                                Value Std.Error DF   t-value p-value
## X(Intercept)                30.47424  2.068067 62 14.735618  0.0000
## Xtrt3control               -10.50819  3.101647 32 -3.387939  0.0019
## Xs(lat_min):trt3cooledFx1   35.46180 18.339913 62  1.933586  0.0577
## Xs(lat_min):trt3controlFx1  12.01946 11.950471 62  1.005773  0.3184
##  Correlation: 
##                            X(Int) Xtrt3c Xs(lt_mn):trt3clF1
## Xtrt3control               -0.667                          
## Xs(lat_min):trt3cooledFx1  -0.009  0.006                   
## Xs(lat_min):trt3controlFx1  0.000  0.016  0.000            
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.00355234 -0.52849216 -0.01912105  0.47693762  4.01379300 
## 
## Number of Observations: 98
## Number of Groups: 
##                    g           g.0 %in% g band %in% g.0 %in% g 
##                    1                    1                   34
# plot(m_cold2, seWithMean = TRUE, shift = coef(m_cold2)[1], residuals = TRUE, rug = FALSE)

m_cold2p <- predict_gam(m_cold2c$gam)
# m_cold2p %>%
#   ggplot(aes(x = lat_min, y = fit)) +
#   geom_smooth_ci(trt3)


# function to compare differences in smooths
# from here:https://fromthebottomoftheheap.net/2017/10/10/difference-splines-i/

smooth_diff <- function(model, newdata, f1, f2, var, alpha = 0.05,
                        unconditional = FALSE) {
    xp <- predict(model, newdata = newdata, type = 'lpmatrix')
    c1 <- grepl(f1, colnames(xp))
    c2 <- grepl(f2, colnames(xp))
    r1 <- newdata[[var]] == f1
    r2 <- newdata[[var]] == f2
    ## difference rows of xp for data from comparison
    X <- xp[r1, ] - xp[r2, ]
    ## zero out cols of X related to splines for other lochs
    X[, ! (c1 | c2)] <- 0
    ## zero out the parametric cols
    X[, !grepl('^s\\(', colnames(xp))] <- 0
    dif <- X %*% coef(model)
    se <- sqrt(rowSums((X %*% vcov(model, unconditional = unconditional)) * X))
    crit <- qt(alpha/2, df.residual(model), lower.tail = FALSE)
    upr <- dif + (crit * se)
    lwr <- dif - (crit * se)
    data.frame(pair = paste(f1, f2, sep = '-'),
               diff = dif,
               se = se,
               upper = upr,
               lower = lwr)
}

# make new simulated data from fit model
pdat <- expand.grid(lat_min = seq(0, 60, by = 1),
                    trt3 = c("cooled", "control"))
xp <- predict(m_cold2c$gam, newdata = pdat, type = 'lpmatrix')

xp2 <- predict(m_cold2c$gam, newdata = pdat)
xp3 <- data.frame(lat_min = pdat$lat_min, trt3 = pdat$trt3, pred = xp2)

comp1 <- smooth_diff(m_cold2c$gam, pdat, "cooled", "control", "trt3")

# the comparison above doesn't account for the parametric estimate of overall group differences
# between cooled and control in the gam. I'm adding it in manually here.

comp1$diff <- comp1$diff + 10.51
comp1$upper <- comp1$upper + 10.51
comp1$lower <- comp1$lower + 10.51

comp1$lat_min <- seq(0, 60, by = 1)

p_exp3 <- ggplot(comp1, mapping = aes(x = lat_min, y = diff)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
  geom_line() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  ylab("Difference in corticosterone") +
  xlab("Sampling latency (minutes)") +
  coord_cartesian(ylim = c(-20, 40)) +
  theme(axis.title = element_text(size = 14), axis.text = element_text(size = 12)) +
  annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "C", size = 6) +
  annotate(geom = "text", x = 30, y = 34, label = "Cooled corticosterone higher", color = "slateblue", size = 4) +
  annotate(geom = "text", x = 30, y = -18, label = "Control corticosterone higher", color = "#CC8400", size = 4)
# read in and combine all raw hobo data
    l <- list.files(here::here("raw_hobo"))
    
    hobo_dat <- data.frame(unit = NA, box = NA, date = NA)
    
    for(i in 1:length(l)){
      lt <- read.csv(here::here("raw_hobo", l[i]))
      lt <- lt[, 2:3]
      colnames(lt) <- c("date", "temp")
      lt$date <- as.POSIXct(lt$date, format = "%m/%d/%Y %H:%M:%S")
      u <- strsplit(l[i], split = "_")[[1]][1]
      b <- strsplit(strsplit(l[i], split = "_")[[1]][2], split = " ")[[1]][[1]]
      lt$unit <- rep(u, nrow(lt))
      lt$box <- rep(b, nrow(lt))
      if(i == 1){hobo_dat <- lt}
      if(i > 1){hobo_dat <- rbind(hobo_dat, lt)}
    }
    
    hobo_dat$yday <- yday(hobo_dat$date)
    hobo_dat$minute <- minute(hobo_dat$date) + 60 * hour(hobo_dat$date)
    hobo_dat$min_2 <- hobo_dat$minute + hobo_dat$yday * 24 * 60
    hobo_dat$ubox <- paste(hobo_dat$unit, hobo_dat$box, sep = "_")

# read in timing of cold snap data
  s <- read.delim(here::here("1_raw_data", "snaps.txt"))
  s$start <- as.POSIXct(s$start, format = "%m/%d/%y %H:%M") 
  s$end <- as.POSIXct(s$end, format = "%m/%d/%y %H:%M")
  s$start2 <- as.POSIXct(s$start2, format = "%m/%d/%y %H:%M") 
  s$end2 <- as.POSIXct(s$end2, format = "%m/%d/%y %H:%M") 

# join hobo and snap data
  hs <- plyr::join(hobo_dat, s, "ubox", "left", "first")
  hs <- subset(hs, hs$experiment == "cold_snap")
  
# loop through hobo and characterize
  hs <- subset(hs, hs$yday >= hs$pend_in & hs$yday <= hs$pend_out)
  
  hs$sn1a <- hs$date - hs$start
  hs$sn1b <- hs$date - hs$end
  
  hs2 <- subset(hs, hs$sn1a > 0 & hs$sn1b < 0)
  hs2 <- subset(hs2, hs2$experiment == "cold_snap")  
  
  hs2$treatment <- gsub("Cold", "Cooled", hs2$treatment)
  p_hobo <- ggplot(data = hs2, mapping = aes(y = temp, x = treatment, color = treatment, fill = treatment)) +
    #geom_jitter(width = 0.15, alpha = 0.9) +
    #geom_boxplot(outlier.shape = NA, alpha = 0.7) + 
    geom_violin(alpha = 0.7) +
    theme_bw() +
    guides(fill = "none", color = "none") +
    ylab("Temperature \u00b0C") +
    xlab("") +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          axis.text.y = element_text(size = 12), axis.title = element_text(size = 14),
          axis.text.x = element_text(size = 12)) +
    geom_point(alpha = 0) +
    stat_summary(fun.data = mean_sd, geom = "pointrange", color = "black", na.rm = TRUE) +
    theme(legend.position = "top") +
    scale_fill_manual(values = c("orange", "slateblue")) +
    scale_color_manual(values = c("orange", "slateblue")) +
    annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "A", size = 6)
  
  
  #p_hobo2 <- ggExtra::ggMarginal(p_hobo, type = 'density', groupFill = TRUE)  
  #ggpubr::ggarrange(p_hobo, p_exp2, p_exp3, nrow = 1)
    
  # hs2$min_5 <- round(hs2$min_2, digits = -1)
  # hs3 <- hs2 %>%
  #   dplyr::group_by(treatment, min_5) %>%
  #   dplyr::summarise(n = n(), temperature = mean(temp)) %>%
  #   pivot_wider(id_cols = min_5, names_from = treatment, values_from = c(temperature, n))
  # 
  # hs3$dif <- hs3$temperature_Cold - hs3$temperature_Control
  # 
  # ggplot(hs3, mapping = aes(x = temperature_Control, y = dif)) +
  #   geom_point(alpha = 0.5) +
  #   geom_smooth(method = "gam")
  
  cowplot::plot_grid(p_hobo, p_exp2, p_exp3, nrow = 1)

2. How and over what timeline does temperature affect natural variation in corticosterone regulation?

For each of basecort, stresscort, and dex I’m fitting a series of sliding window analyses. The full model is identical in every grid except for the choice of temperature time window. Models are general linear mixed models fit with lme4 package like this:

base ~ stage*temperature + (1|ID) + (1|year) induced ~ stage*temperature + base + (1|ID) + (1|year) dex ~ stage*temperature + induced + (1|ID) + (1|year)

2.1 Base Corticosterone

### Basecort heatmap
      
    save_bi <- data.frame(pred = c(lab_cols3, NA), pred_est = NA, ci_lo = NA, ci_hi = NA)  
    save_bp <- data.frame(pred = c(lab_cols3, NA), pred_est = NA, ci_lo = NA, ci_hi = NA)
    for(i in 1:466){
      if(i < 466){
        dbt <- subset(base_af, base_af$stage == "incubation" | base_af$stage == "provisioning")
        dbt$pred <- dbt[, i + 51]
        mbt <- lmer(scale(cort1_correct) ~ scale(pred)*stage + (1|band) + (1|year),
             data = dbt)
        save_bi$pred_est[i] <- fixef(mbt)[2]
        save_bp$pred_est[i] <- fixef(mbt)[2] + fixef(mbt)[4]
        #print(i)
      }
      if(i == 466){
        dbt <- subset(base_af, base_af$stage == "incubation" | base_af$stage == "provisioning")
        dbt$pred <- dbt$t_m_av
        mbt <- lmer(scale(cort1_correct) ~ scale(pred)*stage + (1|band) + (1|year),
             data = dbt)
        save_bi$pred_est[i] <- fixef(mbt)[2]
        save_bp$pred_est[i] <- fixef(mbt)[2] + fixef(mbt)[4]
        #print(i)
      }
      if(i < 466){
        if(str_split(lab_cols3[i], "_")[[1]][2] == str_split(lab_cols3[i], "_")[[1]][3]){
        post <- mvrnorm(n = 1000, mu = fixef(mbt), Sigma = vcov(mbt))
        
        save_bi$ci_lo[i] <- HPDI(post[, 2], prob = 0.95)[1]
        save_bi$ci_hi[i] <- HPDI(post[, 2], prob = 0.95)[2]
        save_bp$ci_lo[i] <- HPDI(post[, 2] + post[, 4], prob = 0.95)[1]
        save_bp$ci_hi[i] <- HPDI(post[, 2] + post[, 4], prob = 0.95)[2]
        }
      }
      if(i == 466){
        post <- mvrnorm(n = 1000, mu = fixef(mbt), Sigma = vcov(mbt))
        
        save_bi$ci_lo[i] <- HPDI(post[, 2], prob = 0.95)[1]
        save_bi$ci_hi[i] <- HPDI(post[, 2], prob = 0.95)[2]
        save_bp$ci_lo[i] <- HPDI(post[, 2] + post[, 4], prob = 0.95)[1]
        save_bp$ci_hi[i] <- HPDI(post[, 2] + post[, 4], prob = 0.95)[2]
      }
      
    }
    
    for(i in 1:(nrow(save_bi) - 1)){
      ssu <- str_split(save_bi$pred[i], "_")[[1]]
      save_bi$start_day[i] <- as.numeric(ssu[2])
      save_bi$period_len[i] <- as.numeric(ssu[3])
      
      ssup <- str_split(save_bp$pred[i], "_")[[1]]
      save_bp$start_day[i] <- as.numeric(ssup[2])
      save_bp$period_len[i] <- as.numeric(ssup[3])
    }
    
      save_bi$start_day[466] <- 0
      save_bi$period_len[466] <- 1
      save_bp$start_day[466] <- 0
      save_bp$period_len[466] <- 1
    
      save_bi7 <- subset(save_bi, save_bi$start_day < 8)
      save_bp7 <- subset(save_bp, save_bp$start_day < 8)
p1 <- ggplot(data = save_bi7, mapping = aes(x = start_day * -1, y = period_len, fill = pred_est)) +
      geom_tile() + 
      scale_fill_gradient2(low = "#542788", high = "#b35806", mid = "#f7f7f7", midpoint = 0,
                           na.value = rgb(211, 211, 211, max = 255, alpha = 0, names = "clear"),
                           limits = c(min(c(save_bi7$pred_est, save_bp7$pred_est)), max(c(save_bi7$pred_est, save_bp7$pred_est)))) +
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 12), axis.title = element_text(size = 12)) +
      theme(legend.position = c(0.85, 0.7)) +
      labs(fill = "Coefficient \n estimate") +
      xlab("Days before capture") +
      ylab("Number of days averaged \n") +
      ggtitle("Incubation: baseline corticosterone") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "A", size = 6) +
      coord_cartesian(ylim = c(0.5, 8), xlim = c(-7.5, 0.5))
    
p2 <- ggplot(data = save_bp7, mapping = aes(x = start_day * -1, y = period_len, fill = pred_est)) +
      geom_tile() + 
      scale_fill_gradient2(low = "#542788", high = "#b35806", mid = "#f7f7f7", midpoint = 0,
                           na.value = rgb(211, 211, 211, max = 255, alpha = 0, names = "clear"),
                           limits = c(min(c(save_bi7$pred_est, save_bp7$pred_est)), max(c(save_bi7$pred_est, save_bp7$pred_est)))) +
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 12), axis.title = element_text(size = 12)) +
      theme(legend.position = c(0.85, 0.7)) +
      labs(fill = "Coefficient \n estimate") +
      xlab("Days before capture") +
      ylab("Number of days averaged \n") +
      ggtitle("Provisioning: baseline corticosterone") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "B", size = 6) +
      coord_cartesian(ylim = c(0.5, 8), xlim = c(-7.5, 0.5))

# do the single line coefficients
    save_bi2 <- subset(save_bi, save_bi$start_day == save_bi$period_len | save_bi$start_day == 0)
    save_bp2 <- subset(save_bp, save_bp$start_day == save_bp$period_len | save_bp$start_day == 0)
    save_bi2 <- subset(save_bi2, save_bi2$start_day < 8)
    save_bp2 <- subset(save_bp2, save_bp2$start_day < 8)
    
p3 <- ggplot(data = save_bi2, mapping = aes(x = start_day * -1, y = pred_est, fill = pred_est)) +
      geom_hline(yintercept = 0, linetype = "dashed") + 
      scale_fill_gradient2(low = "#542788", high = "#b35806", mid = "#f7f7f7", midpoint = 0,
                           na.value = rgb(211, 211, 211, max = 255, alpha = 0, names = "clear"),
                           limits = c(min(c(save_bi7$pred_est, save_bp7$pred_est)), max(c(save_bi7$pred_est, save_bp7$pred_est)))) + 
      geom_line(lwd = 1, color = "gray50") +
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 12), axis.title = element_text(size = 12)) +
      xlab("Days before capture") +
      ylab("Coefficient estimate") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "C", size = 6) +
      #geom_line(data = save_bp2, lwd = 1.5, color = "slateblue") +
      coord_cartesian(ylim = c(-0.3, 0.25), xlim = c(-7.5, 0.5)) +
      #annotate(geom = "text", x = -0.6, y = 0.18, label = "Colder = \n lower cort") +
      #annotate(geom = "text", x = -0.6, y = -0.27, label = "Colder = \n higher cort") +
      #geom_ribbon(mapping = aes(ymin = ci_lo, ymax = ci_hi), fill = "coral3", alpha = 0.4) +
      scale_x_continuous(breaks = seq(-10, 2, 2)) +
      geom_segment(aes(x = start_day * -1, xend = start_day * -1, y = ci_lo, yend = ci_hi), size = 1) +
      geom_point(size = 3, shape = 21, color = "black") +
      guides(fill = "none")

p4 <- ggplot(data = save_bp2, mapping = aes(x = start_day * -1, y = pred_est, fill = pred_est)) +
      geom_hline(yintercept = 0, linetype = "dashed") +
      geom_line(lwd = 1, color = "gray50") + 
      scale_fill_gradient2(low = "#542788", high = "#b35806", mid = "#f7f7f7", midpoint = 0,
                           na.value = rgb(211, 211, 211, max = 255, alpha = 0, names = "clear"),
                           limits = c(min(c(save_bi7$pred_est, save_bp7$pred_est)), max(c(save_bi7$pred_est, save_bp7$pred_est)))) +
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 12), axis.title = element_text(size = 12)) +
      xlab("Days before capture") +
      ylab("Coefficient estimate") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "D", size = 6) +
      coord_cartesian(ylim = c(-0.3, 0.25), xlim = c(-7.5, 0.5)) +
      #annotate(geom = "text", x = -0.6, y = 0.18, label = "Colder = \n lower cort") +
      #annotate(geom = "text", x = -0.6, y = -0.27, label = "Colder = \n higher cort") +
      #geom_ribbon(mapping = aes(ymin = ci_lo, ymax = ci_hi), fill = "coral3", alpha = 0.4) +
      scale_x_continuous(breaks = seq(-10, 2, 2)) +
      geom_segment(aes(x = start_day * -1, xend = start_day * -1, y = ci_lo, yend = ci_hi), size = 1) +
      geom_point(size = 3, color = "black", shape = 21) +
      guides(fill = "none")

ggpubr::ggarrange(p3, p4)
sliding windows for base cort

sliding windows for base cort

## alternative with coefficients in teh same panel
  save_bi2$stage <- "incubation"
  save_bi2$color <- "#009E73"
  save_bi2$start_day2 <- save_bi2$start_day - 0.1
  save_bp2$stage <- "provisioning"
  save_bp2$color <- "#CC79A7"
  save_bp2$start_day2 <- save_bp2$start_day + 0.1
  save_base <- rbind(save_bi2, save_bp2)
  
p_base <-  ggplot(data = save_base, mapping = aes(x = start_day2 * -1, y = pred_est, color = stage)) +
  annotate(geom = "rect", xmin = -0.25, xmax = 0.25, ymin = -0.35, ymax = 0.35, fill = "gray75",
             color = NA, alpha = 0.6) +
  ggtitle("Baseline") +
    geom_hline(yintercept = 0, linetype = "dashed") +
    #geom_line(lwd = 1) +
    theme_bw() +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 14), axis.title = element_text(size = 14)) +
      xlab("Days before capture") +
      ylab("Coefficient estimate") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "A", size = 6) +
      coord_cartesian(ylim = c(-0.4, 0.4), xlim = c(-7.5, 0.5)) +
    scale_x_continuous(breaks = seq(-10, 2, 2)) +
  #scale_y_continuous(breaks = seq(-0.6, 0.6, 0.2)) +
      geom_segment(aes(x = start_day2 * -1, xend = start_day2 * -1, y = ci_lo, yend = ci_hi), size = 1) +
      geom_point(size = 3, fill = "white", shape = 21) +
      guides(fill = "none") +
    theme(legend.position = c(0.55, 0.83), legend.title = element_blank(), 
          legend.box.background = element_rect(fill = "transparent", color = "transparent"),
          legend.text = element_text(size = 12)) +
    scale_color_manual(values = c("#009E5C", "#D38BB3")) #+
    #annotate(geom = "text", x = 0, y = -0.4, label = "*", size = 14)

2.2 Stress-Induced Corticosterone

### stresscort heatmap
      
    save_si <- data.frame(pred = c(lab_cols3, NA), pred_est = NA, ci_lo = NA, ci_hi = NA)  
    save_sp <- data.frame(pred = c(lab_cols3, NA), pred_est = NA, ci_lo = NA, ci_hi = NA)
    for(i in 1:466){
      if(i < 466){
        dst <- subset(stress_af, stress_af$stage == "incubation" | stress_af$stage == "provisioning")
          dst$pred <- dst[, i + 51]
          mst <- lmer(scale(cort2_correct) ~ scale(pred)*stage + scale(cort1_correct) + (1|band) + (1|year),
               data = dst)
          save_si$pred_est[i] <- fixef(mst)[2]
          save_sp$pred_est[i] <- fixef(mst)[2] + fixef(mst)[5]
          #print(i)
      }
      
      if(i == 466){
        dst <- subset(stress_af, stress_af$stage == "incubation" | stress_af$stage == "provisioning")
        dst$pred <- dst$t_m_av
        mst <- lmer(scale(cort2_correct) ~ scale(pred)*stage + scale(cort1_correct) + (1|band) + (1|year),
             data = dst)
        save_si$pred_est[i] <- fixef(mst)[2]
        save_sp$pred_est[i] <- fixef(mst)[2] + fixef(mst)[5]
        #print(i)
      }
          
      if(i < 466){
        if(str_split(lab_cols3[i], "_")[[1]][2] == str_split(lab_cols3[i], "_")[[1]][3]){
        post <- mvrnorm(n = 1000, mu = fixef(mst), Sigma = vcov(mst))
        
        save_si$ci_lo[i] <- HPDI(post[, 2], prob = 0.95)[1]
        save_si$ci_hi[i] <- HPDI(post[, 2], prob = 0.95)[2]
        save_sp$ci_lo[i] <- HPDI(post[, 2] + post[, 5], prob = 0.95)[1]
        save_sp$ci_hi[i] <- HPDI(post[, 2] + post[, 5], prob = 0.95)[2]
        }
      }
      if(i == 466){
        post <- mvrnorm(n = 1000, mu = fixef(mst), Sigma = vcov(mst))
        
        save_si$ci_lo[i] <- HPDI(post[, 2], prob = 0.95)[1]
        save_si$ci_hi[i] <- HPDI(post[, 2], prob = 0.95)[2]
        save_sp$ci_lo[i] <- HPDI(post[, 2] + post[, 5], prob = 0.95)[1]
        save_sp$ci_hi[i] <- HPDI(post[, 2] + post[, 5], prob = 0.95)[2]
      }
        
      
    }
    
    for(i in 1:nrow(save_si)){
      ssu <- str_split(save_si$pred[i], "_")[[1]]
      save_si$start_day[i] <- as.numeric(ssu[2])
      save_si$period_len[i] <- as.numeric(ssu[3])
      
      ssup <- str_split(save_sp$pred[i], "_")[[1]]
      save_sp$start_day[i] <- as.numeric(ssup[2])
      save_sp$period_len[i] <- as.numeric(ssup[3])
    }
    
      save_si$start_day[466] <- 0
      save_si$period_len[466] <- 1
      save_sp$start_day[466] <- 0
      save_sp$period_len[466] <- 1
    
      save_si7 <- subset(save_si, save_si$start_day < 8)
      save_sp7 <- subset(save_sp, save_sp$start_day < 8)
p1s <- ggplot(data = save_si7, mapping = aes(x = start_day * -1, y = period_len, fill = pred_est)) +
      geom_tile() + 
      scale_fill_gradient2(low = "#542788", high = "#b35806", mid = "#f7f7f7", midpoint = 0,
                           na.value = rgb(211, 211, 211, max = 255, alpha = 0, names = "clear"),
                           limits = c(min(c(save_si7$pred_est, save_sp7$pred_est)), max(c(save_si7$pred_est, save_sp7$pred_est)))) +
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 12), axis.title = element_text(size = 12)) +
      theme(legend.position = c(0.85, 0.7)) +
      labs(fill = "Coefficient \n estimate") +
      xlab("Days before capture") +
      ylab("Number of days averaged \n") +
      ggtitle("Incubation: induced corticosterone") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "A", size = 6) +
      coord_cartesian(ylim = c(0.5, 8), xlim = c(-7.5, 0.5))
    
    
p2s <- ggplot(data = save_sp7, mapping = aes(x = start_day * -1, y = period_len, fill = pred_est)) +
      geom_tile() + 
      scale_fill_gradient2(low = "#542788", high = "#b35806", mid = "#f7f7f7", midpoint = 0,
                           na.value = rgb(211, 211, 211, max = 255, alpha = 0, names = "clear"),
                           limits = c(min(c(save_si7$pred_est, save_sp7$pred_est)), max(c(save_si7$pred_est, save_sp7$pred_est)))) +
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 12), axis.title = element_text(size = 12)) +
      theme(legend.position = c(0.85, 0.7)) +
      labs(fill = "Coefficient \n estimate") +
      xlab("Days before capture") +
      ylab("Number of days averaged \n") +
      ggtitle("Provisioning: induced corticosterone") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "B", size = 6) +
      coord_cartesian(ylim = c(0.5, 8), xlim = c(-7.5, 0.5))

# do the single line coefficients
    save_si2 <- subset(save_si, save_si$start_day == save_si$period_len | save_si$start_day == 0)
    save_sp2 <- subset(save_sp, save_sp$start_day == save_sp$period_len | save_sp$start_day == 0)
    save_si2 <- subset(save_si2, save_si2$start_day < 8)
    save_sp2 <- subset(save_sp2, save_sp2$start_day < 8)
    
p3s <- ggplot(data = save_si2, mapping = aes(x = start_day * -1, y = pred_est, fill = pred_est)) +
      geom_hline(yintercept = 0, linetype = "dashed") +
      geom_line(lwd = 1, color = "gray50") + 
      scale_fill_gradient2(low = "#542788", high = "#b35806", mid = "#f7f7f7", midpoint = 0,
                           na.value = rgb(211, 211, 211, max = 255, alpha = 0, names = "clear"),
                           limits = c(min(c(save_si7$pred_est, save_sp7$pred_est)), max(c(save_si7$pred_est, save_sp7$pred_est)))) +
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 12), axis.title = element_text(size = 12)) +
      xlab("Days before capture") +
      ylab("Coefficient estimate") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "C", size = 6) +
      #geom_line(data = save_sp2, lwd = 1.5, color = "slateblue") +
      coord_cartesian(xlim = c(-7.5, 0.5), ylim = c(-0.4, 0.15)) +
      #annotate(geom = "text", x = 2, y = 0.15, label = "Warmer = \n higher cort") +
      #annotate(geom = "text", x = 2, y = -0.08, label = "Warmer = \n lower cort") +
      #geom_ribbon(mapping = aes(ymin = ci_lo, ymax = ci_hi), fill = "coral3", alpha = 0.4) +
      scale_x_continuous(breaks = seq(-10, 2, 2)) +
      geom_segment(aes(x = start_day * -1, xend = start_day * -1, y = ci_lo, yend = ci_hi), size = 1) +
      geom_point(size = 3, color = "black", shape = 21) +
      guides(fill = "none")

p4s <- ggplot(data = save_sp2, mapping = aes(x = start_day * -1, y = pred_est, fill = pred_est)) +
      geom_hline(yintercept = 0, linetype = "dashed") +
      geom_line(lwd = 1, color = "gray50") + 
      scale_fill_gradient2(low = "#542788", high = "#b35806", mid = "#f7f7f7", midpoint = 0,
                           na.value = rgb(211, 211, 211, max = 255, alpha = 0, names = "clear"),
                           limits = c(min(c(save_si7$pred_est, save_sp7$pred_est)), max(c(save_si7$pred_est, save_sp7$pred_est)))) +
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 12), axis.title = element_text(size = 12)) +
      xlab("Days before capture") +
      ylab("Coefficient estimate") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "D", size = 6) +
      coord_cartesian(xlim = c(-7.5, 0.5), ylim = c(-0.4, 0.15)) +
      #annotate(geom = "text", x = 2, y = 0.15, label = "Warmer = \n higher cort") +
      #annotate(geom = "text", x = 2, y = -0.08, label = "Warmer = \n lower cort") +
      #geom_ribbon(mapping = aes(ymin = ci_lo, ymax = ci_hi), fill = "coral3", alpha = 0.4) +
      scale_x_continuous(breaks = seq(-10, 2, 2)) +
      geom_segment(aes(x = start_day * -1, xend = start_day * -1, y = ci_lo, yend = ci_hi), size = 1) +
      geom_point(size = 3, color = "black", shape = 21) +
      guides(fill = "none")

ggpubr::ggarrange(p3s, p4s)
sliding windows for stress

sliding windows for stress

## alternative with coefficients in teh same panel
  save_si2$stage <- "incubation"
  save_si2$color <- "#009E73"
  save_si2$start_day2 <- save_si2$start_day - 0.1
  save_sp2$stage <- "provisioning"
  save_sp2$color <- "#CC79A7"
  save_sp2$start_day2 <- save_sp2$start_day + 0.1
  save_stress <- rbind(save_si2, save_sp2)
  
p_stress <-  ggplot(data = save_stress, mapping = aes(x = start_day2 * -1, y = pred_est, color = stage)) +
    #geom_rect(mapping = aes(xmin = -0.5, xmax = 0.5, ymin = -0.35, ymax = 0.35, fill = t), 
    #          fill = "goldenrod", alpha = 0.2, color = NA) +
  ggtitle("Stress-induced") +
    annotate(geom = "rect", xmin = -3.25, xmax = -2.75, ymin = -0.35, ymax = 0.35, fill = "gray75",
             color = NA, alpha = 0.6) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    #geom_line(lwd = 1) +
    theme_bw() +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 14), axis.title = element_text(size = 14)) +
      xlab("Days before capture") +
      ylab("Coefficient estimate") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "B", size = 6) +
      coord_cartesian(xlim = c(-7.5, 0.5), ylim = c(-0.4, 0.4)) +
    scale_x_continuous(breaks = seq(-10, 2, 2)) +
  #scale_y_continuous(breaks = seq(-0.6, 0.6, 0.2)) +
      geom_segment(aes(x = start_day2 * -1, xend = start_day2 * -1, y = ci_lo, yend = ci_hi), size = 1) +
      geom_point(size = 3, fill = "white", shape = 21) +
      guides(fill = "none", color = "none") +
    #theme(legend.position = c(0.8, 0.85), legend.title = element_blank()) +
    scale_color_manual(values = c("#009E5C", "#D38BB3")) 
    #annotate(geom = "text", x = -3, y = -0.4, label = "*", size = 14)
#p_stress

2.3 Post-Dex Corticosterone

### stresscort heatmap
      
    save_di <- data.frame(pred = c(lab_cols3, NA), pred_est = NA, ci_lo = NA, ci_hi = NA)  
    save_dp <- data.frame(pred = c(lab_cols3, NA), pred_est = NA, ci_lo = NA, ci_hi = NA)
    for(i in 1:466){
      if(i < 466){
        ddt <- subset(dex_af, dex_af$stage == "incubation" | dex_af$stage == "provisioning")
        ddt$pred <- ddt[, i + 51]
        mdt <- lmer(scale(cort3_correct) ~ scale(pred)*stage + scale(cort2_correct) + (1|band) + (1|year),
             data = ddt)
        save_di$pred_est[i] <- fixef(mdt)[2]
        save_dp$pred_est[i] <- fixef(mdt)[2] + fixef(mdt)[5]
        #print(i)
      }
      if(i == 466){
        ddt <- subset(dex_af, dex_af$stage == "incubation" | dex_af$stage == "provisioning")
        ddt$pred <- ddt$t_m_av
        mdt <- lmer(scale(cort3_correct) ~ scale(pred)*stage + scale(cort2_correct) + (1|band) + (1|year),
             data = ddt)
        save_di$pred_est[i] <- fixef(mdt)[2]
        save_dp$pred_est[i] <- fixef(mdt)[2] + fixef(mdt)[5]
        #print(i)
      }
      
      if(i < 466){
        if(str_split(lab_cols3[i], "_")[[1]][2] == str_split(lab_cols3[i], "_")[[1]][3]){
        post <- mvrnorm(n = 1000, mu = fixef(mdt), Sigma = vcov(mdt))
        
        save_di$ci_lo[i] <- HPDI(post[, 2], prob = 0.95)[1]
        save_di$ci_hi[i] <- HPDI(post[, 2], prob = 0.95)[2]
        save_dp$ci_lo[i] <- HPDI(post[, 2] + post[, 5], prob = 0.95)[1]
        save_dp$ci_hi[i] <- HPDI(post[, 2] + post[, 5], prob = 0.95)[2]
        }
      }
      if(i == 466){
        post <- mvrnorm(n = 1000, mu = fixef(mdt), Sigma = vcov(mdt))
        
        save_di$ci_lo[i] <- HPDI(post[, 2], prob = 0.95)[1]
        save_di$ci_hi[i] <- HPDI(post[, 2], prob = 0.95)[2]
        save_dp$ci_lo[i] <- HPDI(post[, 2] + post[, 5], prob = 0.95)[1]
        save_dp$ci_hi[i] <- HPDI(post[, 2] + post[, 5], prob = 0.95)[2]
      }  
        
        
      
    }
    
    for(i in 1:nrow(save_di)){
      ssu <- str_split(save_di$pred[i], "_")[[1]]
      save_di$start_day[i] <- as.numeric(ssu[2])
      save_di$period_len[i] <- as.numeric(ssu[3])
      
      ssup <- str_split(save_dp$pred[i], "_")[[1]]
      save_dp$start_day[i] <- as.numeric(ssup[2])
      save_dp$period_len[i] <- as.numeric(ssup[3])
    }
    
      save_di$start_day[466] <- 0
      save_di$period_len[466] <- 1
      save_dp$start_day[466] <- 0
      save_dp$period_len[466] <- 1
    
      
      save_di21 <- subset(save_di, save_di$start_day < 22)
      save_dp21 <- subset(save_dp, save_dp$start_day < 22)
p1d <- ggplot(data = save_di21, mapping = aes(x = start_day * -1, y = period_len, fill = pred_est)) +
      geom_tile() + 
      scale_fill_gradient2(low = "#542788", high = "#b35806", mid = "#f7f7f7", midpoint = 0,
                           na.value = rgb(211, 211, 211, max = 255, alpha = 0, names = "clear"),
                           limits = c(min(c(save_di21$pred_est, save_dp21$pred_est)), max(c(save_di21$pred_est, save_dp21$pred_est)))) +
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 12), axis.title = element_text(size = 12)) +
      theme(legend.position = c(0.85, 0.7)) +
      labs(fill = "Coefficient \n estimate") +
      xlab("Days before capture") +
      ylab("Number of days averaged \n") +
      ggtitle("Incubation: post-dex corticosterone") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "A", size = 6) +
      coord_cartesian(ylim = c(0.5, 23))
    
p2d <- ggplot(data = save_dp21, mapping = aes(x = start_day * -1, y = period_len, fill = pred_est)) +
      geom_tile() + 
      scale_fill_gradient2(low = "#542788", high = "#b35806", mid = "#f7f7f7", midpoint = 0,
                           na.value = rgb(211, 211, 211, max = 255, alpha = 0, names = "clear"),
                           limits = c(min(c(save_di21$pred_est, save_dp21$pred_est)), max(c(save_di21$pred_est, save_dp21$pred_est)))) +
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 12), axis.title = element_text(size = 12)) +
      theme(legend.position = c(0.85, 0.7)) +
      labs(fill = "Coefficient \n estimate") +
      xlab("Days before capture") +
      ylab("Number of days averaged \n") +
      ggtitle("Provisioning: post-dex corticosterone") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "B", size = 6) +
      coord_cartesian(ylim = c(0.5, 23))

# do the single line coefficients
    save_di2 <- subset(save_di, save_di$start_day == save_di$period_len | save_di$start_day == 0)
    save_dp2 <- subset(save_dp, save_dp$start_day == save_dp$period_len | save_dp$start_day == 0)
    save_di2 <- subset(save_di2, save_di2$start_day < 22)
    save_dp2 <- subset(save_dp2, save_dp2$start_day < 22)
    
p3d <- ggplot(data = save_di2, mapping = aes(x = start_day * -1, y = pred_est, fill = pred_est)) +
      geom_hline(yintercept = 0, linetype = "dashed") +
      geom_line(lwd = 1, color = "gray50") + 
      scale_fill_gradient2(low = "#542788", high = "#b35806", mid = "#f7f7f7", midpoint = 0,
                           na.value = rgb(211, 211, 211, max = 255, alpha = 0, names = "clear"),
                           limits = c(min(c(save_di21$pred_est, save_dp21$pred_est)), max(c(save_di21$pred_est, save_dp21$pred_est)))) +
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 12), axis.title = element_text(size = 12)) +
      xlab("Days before capture") +
      ylab("Coefficient estimate") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "C", size = 6) +
      #geom_line(data = save_dp2, lwd = 1.5, color = "slateblue") +
      coord_cartesian(ylim = c(-0.25, 0.35)) +
      #annotate(geom = "text", x = 2, y = 0.15, label = "Warmer = \n higher cort") +
      #annotate(geom = "text", x = 2, y = -0.08, label = "Warmer = \n lower cort") +
      #geom_ribbon(mapping = aes(ymin = ci_lo, ymax = ci_hi), fill = "coral3", alpha = 0.4) +
      scale_x_continuous(breaks = seq(-30, 5, 5)) +
      geom_segment(aes(x = start_day * -1, xend = start_day * -1, y = ci_lo, yend = ci_hi), size = 1) +
      geom_point(size = 3, color = "black", shape = 21) +
      guides(fill = "none")

p4d <- ggplot(data = save_dp2, mapping = aes(x = start_day * -1, y = pred_est, fill = pred_est)) +
      geom_hline(yintercept = 0, linetype = "dashed") +
      geom_line(lwd = 1, color = "gray50") + 
      scale_fill_gradient2(low = "#542788", high = "#b35806", mid = "#f7f7f7", midpoint = 0,
                           na.value = rgb(211, 211, 211, max = 255, alpha = 0, names = "clear"),
                           limits = c(min(c(save_di21$pred_est, save_dp21$pred_est)), max(c(save_di21$pred_est, save_dp21$pred_est)))) +
      theme_bw() +
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 12), axis.title = element_text(size = 12)) +
      xlab("Days before capture") +
      ylab("Coefficient estimate") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "D", size = 6) +
      coord_cartesian(ylim = c(-0.25, 0.35)) +
      #annotate(geom = "text", x = 2, y = 0.15, label = "Warmer = \n higher cort") +
      #annotate(geom = "text", x = 2, y = -0.08, label = "Warmer = \n lower cort") +
      #geom_ribbon(mapping = aes(ymin = ci_lo, ymax = ci_hi), fill = "coral3", alpha = 0.4) +
      scale_x_continuous(breaks = seq(-30, 5, 5)) +
      geom_segment(aes(x = start_day * -1, xend = start_day * -1, y = ci_lo, yend = ci_hi), size = 1) +
      geom_point(size = 3, color = "black", shape = 21) +
      guides(fill = "none")

ggpubr::ggarrange(p3d, p4d)
sliding windows for dex

sliding windows for dex

## alternative with coefficients in teh same panel
  save_di2$stage <- "incubation"
  save_di2$color <- "#009E73"
  save_di2$start_day2 <- save_di2$start_day - 0.2
  save_dp2$stage <- "provisioning"
  save_dp2$color <- "#CC79A7"
  save_dp2$start_day2 <- save_dp2$start_day + 0.2
  save_dex <- rbind(save_di2, save_dp2)
  
p_dex <-  ggplot(data = save_dex, mapping = aes(x = start_day2 * -1, y = pred_est, color = stage)) +
  annotate(geom = "rect", xmin = -0.3, xmax = 0.3, ymin = -0.35, ymax = 0.35, fill = "gray75",
             color = NA, alpha = 0.6) +
  ggtitle("Post-dexamethasone") +
    geom_hline(yintercept = 0, linetype = "dashed") +
    #geom_line(lwd = 1) +
    theme_bw() +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
            axis.text = element_text(size = 14), axis.title = element_text(size = 14)) +
      xlab("Days before capture") +
      ylab("Coefficient estimate") +
      annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "C", size = 6) +
      coord_cartesian(ylim = c(-0.4, 0.4)) +
    scale_x_continuous(breaks = seq(-30, 5, 5)) +
  #scale_y_continuous(breaks = seq(-0.6, 0.6, 0.2)) +
      geom_segment(aes(x = start_day2 * -1, xend = start_day2 * -1, y = ci_lo, yend = ci_hi), size = 1) +
      geom_point(size = 3, fill = "white", shape = 21) +
      guides(fill = "none", color = "none") +
    #theme(legend.position = c(0.8, 0.85), legend.title = element_blank()) +
    scale_color_manual(values = c("#009E5C", "#D38BB3")) #+
    #annotate(geom = "text", x = 0, y = -0.4, label = "*", size = 14)
#p_dex

2.4 Combined plot

This is the same as the panels above but put together into the version included with the paper.

3. Is variation in corticosterone best predicted by absolute temperature or a prediction error?

Using AIC to compare models that have either absolute temperature or prediction error (real vs. expected temp based on date). See text for description. For each model set I am printing the AIC table and the details for best supported model(s).

3.1 Base Corticosterone

##               df   logLik    AICc  delta weight
## b_temp_stg     8 -1811.60 3639.32  0.000  0.900
## b_temp         6 -1816.54 3645.13  5.820  0.049
## b_offs_stg     8 -1814.79 3645.68  6.365  0.037
## b_temp_stg_db 12 -1812.50 3649.22  9.909  0.006
## b_null         6 -1818.98 3650.02 10.709  0.004
## b_offs         6 -1819.40 3650.86 11.542  0.003
## b_temp_db      8 -1818.66 3653.42 14.107  0.001
  cort 1 correct
Predictors Estimates CI p
(Intercept) 0.75 0.45 – 1.06 <0.001
temperature -0.14 -0.20 – -0.09 <0.001
stage [provisioning] 0.27 0.13 – 0.40 <0.001
cap_doy -0.07 -0.14 – -0.00 0.044
temperature * stage
[provisioning]
0.08 -0.02 – 0.18 0.123
Random Effects
σ2 0.74
τ00 band 0.07
τ00 year 0.18
ICC 0.25
N year 8
N band 689
Observations 1365
Marginal R2 / Conditional R2 0.026 / 0.271

3.2 Stress Cort

##               df   logLik    AICc  delta weight
## s_temp_stg     9 -1033.79 2085.76  0.000  0.951
## s_temp         7 -1039.39 2092.89  7.136  0.027
## s_offs_stg     9 -1037.58 2093.34  7.580  0.021
## s_offs         7 -1042.93 2099.97 14.211  0.001
## s_temp_db      9 -1049.10 2116.38 30.622  0.000
## s_temp_stg_db 13 -1045.11 2116.59 30.832  0.000
## s_null         7 -1057.88 2129.87 44.117  0.000
  cort 2 correct
Predictors Estimates CI p
(Intercept) 3.31 3.17 – 3.46 <0.001
temperature -0.19 -0.24 – -0.14 <0.001
stage [provisioning] -0.20 -0.31 – -0.09 <0.001
cap_doy -0.04 -0.10 – 0.02 0.169
cort1_correct 0.11 0.06 – 0.15 <0.001
temperature * stage
[provisioning]
0.12 0.04 – 0.20 0.003
Random Effects
σ2 0.24
τ00 band 0.27
τ00 year 0.04
ICC 0.56
N year 8
N band 589
Observations 1011
Marginal R2 / Conditional R2 0.097 / 0.606

3.3 Dex Cort

##               df  logLik   AICc  delta weight
## d_null         7 -321.24 656.71  0.000  0.865
## d_temp         7 -323.68 661.58  4.874  0.076
## d_offs         7 -323.96 662.14  5.433  0.057
## d_temp_stg     9 -325.45 669.26 12.551  0.002
## d_offs_stg     9 -326.05 670.45 13.745  0.001
## d_temp_db      9 -327.41 673.19 16.481  0.000
## d_temp_stg_db 13 -332.25 691.24 34.531  0.000
  cort 3 correct
Predictors Estimates CI p
(Intercept) 2.69 2.46 – 2.92 <0.001
stage [provisioning] -0.14 -0.27 – -0.00 0.049
cap_doy 0.03 -0.02 – 0.09 0.196
cort2_correct 0.18 0.14 – 0.22 <0.001
Random Effects
σ2 0.19
τ00 band 0.00
τ00 year 0.09
ICC 0.33
N year 7
N band 346
Observations 512
Marginal R2 / Conditional R2 0.114 / 0.409

4. Plots from best models

4.1 Baseline plot

For these plots I am taking the best model from the section above and calculating the mean and confidence interval from the table of coefficients separately for incubation and provisioning. I do this by taking 10000 samples from the posterior distribution and calculating CIs at each point then plotting. Models are fit on log scale and all CIs/means are also calculated on that scale, but for plotting I back-transform to the measurement scale to make it easier to interpret. Full details on the models themselves are given in the tables above.

postb <- mvrnorm(n = 10000, mu = fixef(b_temp_stg), Sigma = vcov(b_temp_stg))
rtb <- seq(min(base_af_comp$temperature), max(base_af_comp$temperature), length.out = 100)

b_mu_inc <- exp(sapply(rtb, function(z)mean(postb[, 1] + postb[, 2]*z)))
b_ci_inc <- exp(sapply(rtb, function(z)HPDI(postb[, 1] + postb[, 2]*z)))

b_mu_pro <- exp(sapply(rtb, function(z)mean(postb[, 1] + postb[, 2]*z + postb[, 3] + postb[, 5]*z)))
b_ci_pro <- exp(sapply(rtb, function(z)HPDI(postb[, 1] + postb[, 2]*z + postb[, 3] + postb[, 5]*z)))

rtb2 <- (rtb * sd(na.omit(base_af$t_m_av))) + mean(na.omit(base_af$t_m_av))

b_mu <- data.frame(temp = rtb2, mu_inc = b_mu_inc, mu_pro = b_mu_pro)
b_ci <- data.frame(temp = rtb2, cil_inc = b_ci_inc[1,], cih_inc = b_ci_inc[2,], cil_pro = b_ci_pro[1,], cih_pro = b_ci_pro[2,])


p2_base <- ggplot(data = base_af_comp, mapping = aes(x = (temperature * sd(na.omit(base_af$t_m_av))) + mean(na.omit(base_af$t_m_av)), 
                                            y = exp(cort1_correct), fill = stage)) +
  geom_point(alpha = 0.5, shape = 21) +
  scale_fill_manual(values = c("#009E5C", "#D38BB3")) +
  geom_line(data = b_mu, mapping = aes(x = temp, y = mu_inc), inherit.aes = FALSE, color = "#009E5C", size = 1.5, alpha = 0.9) +
  geom_line(data = b_mu, mapping = aes(x = temp, y = mu_pro), inherit.aes = FALSE, color = "#D38BB3", size = 1.5, alpha = 0.9) +
  geom_ribbon(data = b_ci, mapping = aes(x = temp, ymin = cil_inc, ymax = cih_inc), inherit.aes = FALSE, fill = "#009E5C", alpha = 0.3) +
  geom_ribbon(data = b_ci, mapping = aes(x = temp, ymin = cil_pro, ymax = cih_pro), inherit.aes = FALSE, fill = "#D38BB3", alpha = 0.3) +
  theme_bw() +
  annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "A", size = 6) +
  coord_cartesian(ylim = c(0, 15)) +
  theme(legend.position = c(0.81, 0.91), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), 
        axis.title = element_text(size = 16), legend.title = element_blank(), axis.text = element_text(size = 14),
        legend.text = element_text(size = 14)) +
  guides(fill = guide_legend(override.aes = list(size = 4))) +
  ylab("Baseline corticosterone (ng/ml)") +
  xlab("Temperature \u00b0C") 
p2_base

Note that in this plot I’ve clipped the y axis to 20 to make the trend visible at all. This means 34 data points (out of 1428) are not plotted.

4.2 Stress induced plot

Exact same approach as described for baseline plot above.

posts <- mvrnorm(n = 10000, mu = fixef(s_temp_stg), Sigma = vcov(s_temp_stg))
rt <- seq(min(stress_af_comp$temperature), max(stress_af_comp$temperature), length.out = 100)

s_mu_inc <- exp(sapply(rt, function(z)mean(posts[, 1] + posts[, 2]*z)))
s_ci_inc <- exp(sapply(rt, function(z)HPDI(posts[, 1] + posts[, 2]*z)))

s_mu_pro <- exp(sapply(rt, function(z)mean(posts[, 1] + posts[, 2]*z + posts[, 3] + posts[, 6]*z)))
s_ci_pro <- exp(sapply(rt, function(z)HPDI(posts[, 1] + posts[, 2]*z + posts[, 3] + posts[, 6]*z)))

rt2 <- (rt * sd(na.omit(stress_af$tb_3_3))) + mean(na.omit(stress_af$tb_3_3))

s_mu <- data.frame(temp = rt2, mu_inc = s_mu_inc, mu_pro = s_mu_pro)
s_ci <- data.frame(temp = rt2, cil_inc = s_ci_inc[1,], cih_inc = s_ci_inc[2,], cil_pro = s_ci_pro[1,], cih_pro = s_ci_pro[2,])

#pred <- predict(ss_temp_stg)

p2_stress <- ggplot(data = stress_af_comp, mapping = aes(x = (temperature * sd(na.omit(stress_af$tb_3_3))) + mean(na.omit(stress_af$tb_3_3)), 
                                            y = exp(cort2_correct), fill = stage)) +
  geom_point(alpha = 0.5, shape = 21) +
  scale_fill_manual(values = c("#009E5C", "#D38BB3")) +
  geom_line(data = s_mu, mapping = aes(x = temp, y = mu_inc), inherit.aes = FALSE, color = "#009E5C", size = 1.5, alpha = 0.9) +
  geom_line(data = s_mu, mapping = aes(x = temp, y = mu_pro), inherit.aes = FALSE, color = "#D38BB3", size = 1.5, alpha = 0.9) +
  geom_ribbon(data = s_ci, mapping = aes(x = temp, ymin = cil_inc, ymax = cih_inc), inherit.aes = FALSE, fill = "#009E5C", alpha = 0.3) +
  geom_ribbon(data = s_ci, mapping = aes(x = temp, ymin = cil_pro, ymax = cih_pro), inherit.aes = FALSE, fill = "#D38BB3", alpha = 0.3) +
  theme_bw() +
  coord_cartesian(ylim = c(0, 120)) +
  annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "B", size = 6) +
  theme(legend.position = c(0.81, 0.91), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), 
        axis.title = element_text(size = 16), legend.title = element_blank(), axis.text = element_text(size = 14),
        legend.text = element_text(size = 14)) +
  guides(fill = guide_legend(override.aes = list(size = 4))) +
  ylab("Stress-induced corticosterone (ng/ml)") +
  xlab("Temperature \u00b0C") 
p2_stress

The two panels for base and stress-induced combined into one figure for the manuscript.

4.3 Dex plot

The best model is the null with no temperature effect so not plot is produced.